THE  WIENER-ASKEY  POLYNOMIAL  CHAOS  FOR  STOCHASTIC 
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Abstract.  We  present  a  new  method  for  solving  stochastic  differential  equations  based  on 
Galerkin  projections  and  extensions  of  Wiener’s  polynomial  chaos.  Specifically,  we  represent  the 
stochastic  processes  with  an  optimum  trial  basis  from  the  Askey  family  of  orthogonal  polynomials 
that  reduces  the  dimensionality  of  the  system  and  leads  to  exponential  convergence  of  the  error. 
Several  continuous  and  discrete  processes  are  treated,  and  numerical  examples  show  substantial 
speed-up  compared  to  Monte-Carlo  simulations  for  low  dimensional  stochastic  inputs. 
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1.  Introduction.  Wiener  first  defined  ‘Homogeneous  Chaos’  as  the  span  of  Her- 
mite  polynomial  functionals  of  a  Gaussian  process  [19];  Polynomial  Chaos  is  defined  as 
the  member  of  that  set.  According  to  the  Cameron-Martin  theorem  [3],  the  Fourier- 
Hermite  series  converge  to  any  functional  in  the  sense.  In  the  context  of 
stochastic  processes,  this  implies  that  the  homogeneous  chaos  expansion  converges 
to  any  processes  with  finite  second-order  moments.  Therefore,  such  an  expansion 
provides  a  means  of  representing  a  stochastic  process  with  Hermite  orthogonal  poly¬ 
nomials.  Other  names  such  as  ‘Wiener  chaos’,  ‘Wiener-Hermite  chaos’,  etc.,  have  also 
been  used  in  the  literature.  In  this  paper,  we  will  use  the  term  Hermite-Chaos. 

While  Hermite-Chaos  is  useful  in  the  analysis  of  stochastic  processes,  efforts  have 
also  been  made  to  apply  it  to  model  uncertainty  in  physical  applications.  In  this 
case,  the  continuous  integral  form  of  the  Hermite-Chaos  is  written  in  the  discrete 
form  of  infinite  summation  which  is  further  truncated.  Ghanem  &  Spanos  [9]  com¬ 
bined  the  Hermite-Chaos  expansion  with  finite  element  method  to  model  uncertainty 
encountered  in  various  problems  of  solid  mechanics,  e.g.,  [7],  [8],  [9],  etc.  In  [20],  the 
polynomial  chaos  was  applied  to  modeling  uncertainty  in  fluid  dynamics  applications. 
The  algorithm  was  implemented  in  the  context  of  spectral/ft.p  element  method  and 
various  benchmark  tests  were  conducted  to  demonstrate  convergence  in  prototype 
flows. 

Although  for  any  arbitrary  random  process  with  finite  second-order  moments,  the 
Hermite-Chaos  expansion  converges  in  accord  with  Cameron-Martin  theorem  [3],  it 
has  been  demonstrated  that  the  convergence  rate  is  optimal  for  Gaussian  processes, 
in  fact  the  rate  is  exponential  [15].  This  can  be  understood  from  the  fact  that  the 
weighting  function  of  Hermite  polynomials  is  the  same  as  the  probability  density  func¬ 
tion  of  the  Gaussian  random  variables.  For  other  types  of  processes  the  convergence 
rate  may  be  substantially  slower.  In  this  case,  other  types  of  orthogonal  polynomials, 
instead  of  Hermite  polynomials,  could  be  used  to  construct  the  chaos  expansion.  In 
an  early  work  by  Ogura  [16],  a  chaos  expansion  based  on  Charlier  polynomials  was 
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proposed  to  represent  the  Poisson  processes,  following  the  theory  of  ‘discrete  chaos’ 
by  Wiener  [19]. 

An  important  class  of  orthogonal  polynomials  are  the  members  of  the  so-called 
Askey-scheme  of  polynomials  [1].  This  scheme  classifies  the  hypergeometric  orthog¬ 
onal  polynomials  that  satisfy  some  type  of  differential  or  difference  equation  and 
indicates  the  limit  relations  between  them.  Hermite  polynomials  are  a  subset  of  the 
Askey-scheme.  Each  subset  of  the  orthogonal  polynomials  in  the  Askey-scheme  has 
different  weighting  function  in  their  orthogonality  relationship.  It  has  been  realized 
that  some  of  these  weighting  functions  are  identical  to  the  probability  function  of 
certain  random  distributions.  For  example, 

•  Hermite  polynomials  are  associated  with  the  Gaussian  distribution. 

•  Laguerre  polynomials  with  the  Gamma  distribution. 

•  Jacobi  polynomials  with  the  Beta  distribution. 

•  Gharlier  polynomials  with  the  Poisson  distribution. 

•  Meixner  polynomials  with  the  negative  Binomial  distribution. 

•  Krawtchouk  polynomials  with  the  Binomial  distribution,  and 

•  Hahn  polynomials  with  the  Hypergeometric  distribution. 

This  finding  opens  the  possibility  of  representing  stochastic  processes  with  different 
orthogonal  polynomials  according  to  the  property  of  the  processes. 

The  close  connection  between  stochastic  processes  and  orthogonal  polynomials  has 
long  been  recognized.  Despite  of  the  role  of  Hermite  polynomials  in  the  integration 
theory  of  Brownian  motion  ([19]  and  [11]),  many  birth-and-death  models  were  related 
to  specific  orthogonal  polynomials.  The  so-called  Karlin-McGregor  representation 
of  the  transition  probabilities  of  a  birth-and-death  process  is  in  terms  of  orthogonal 
polynomials  [12].  In  [16]  and  [5],  the  integral  relation  between  the  Poisson  process  and 
the  Gharlier  polynomials  was  found.  In  [17]  the  role  of  the  orthogonal  polynomials 
from  Askey-scheme  in  the  theory  of  Markov  processes  was  studied  and  the  connection 
between  the  Krawtchouk  polynomials  and  the  binomial  process  was  established. 

In  this  paper,  we  extend  the  work  by  Ghanem  &  Spanos  for  Hermite-Ghaos  expan¬ 
sion  [9]  and  Ogura  for  Gharlier-Ghaos  expansion  [16].  We  propose  an  Askey-scheme 
based  polynomial  chaos  expansion  for  stochastic  processes  which  includes  all  the  or¬ 
thogonal  polynomials  in  the  above  list.  We  demonstrate  numerically  the  optimal 
(exponential)  convergence  rate  of  each  Wiener- Askey  polynomial  chaos  expansion  for 
their  corresponding  stochastic  processes  by  solving  a  stochastic  ordinary  differential 
equation,  for  which  the  exact  solutions  can  be  obtained.  It  is  also  shown  that  if 
for  a  certain  process  the  optimal  Wiener-Askey  polynomial  chaos  expansion  is  not 
employed,  the  solution  also  converges  but  the  rate  is  clearly  slower.  This  approach 
will  provide  a  guideline  for  representing  stochastic  processes  in  physical  applications 
properly. 

In  practical  applications,  one  often  does  not  know  the  analytical  form  of  the  dis¬ 
tribution  of  the  process,  or  if  known,  it  may  not  be  one  of  the  basic  distributions,  e.g. 
Gaussian,  Poisson,  etc..  In  this  case,  one  can  choose  a  set  of  Wiener-Askey  polyno¬ 
mial  chaos  expansion  and  conduct  a  numerical  projection  procedure  to  represent  the 
process.  This  issue  will  be  addressed  in  the  present  paper  as  well. 

This  paper  is  organized  as  follows:  In  the  next  section  we  review  the  theory  of  the 
Askey-scheme  of  hypergeometric  orthogonal  polynomials,  and  in  section  3  we  review 
the  theory  of  the  original  Wiener’s  polynomial  chaos.  In  section  4  we  present  the 
framework  of  Wiener-Askey  polynomial  chaos  expansion  for  stochastic  processes.  In 
section  5  we  present  numerical  solutions  of  a  stochastic  ordinary  differential  equation 
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with  different  Wiener-Askey  chaos  expansions.  The  choice  of  the  particular  Wiener- 
Askey  chaos  is  based  on  the  distribution  of  the  random  input,  and  we  demonstrate  the 
exponential  convergence  rate  with  the  appropriately  chosen  Wiener-Askey  basis.  In 
section  6,  we  address  the  issue  of  representing  an  arbitrary  random  distribution  and 
we  show  that  although  the  Wiener-Askey  polynomial  chaos  converges  in  general,  the 
exponential  convergence  is  not  realized  if  the  optimal  type  of  Wiener-Askey  chaos  is 
not  chosen.  We  conclude  the  paper  with  a  discussion  on  possible  extensions  and  appli¬ 
cations  to  more  complicated  problems.  An  appendix  of  the  definitions  and  properties 
of  the  orthogonal  polynomials  discussed  in  this  paper  is  included  for  completeness. 

2.  The  Askey-scheme  of  Hypergeometric  Orthogonal  Polynomials.  The 

theory  of  orthogonal  polynomials  is  relatively  mature  and  several  books  have  been 
devoted  to  their  study  (e.g.,  [18],  [2],  [4],  etc.).  However,  more  recent  work  has 
shown  that  an  important  class  of  orthogonal  polynomials  belong  to  the  Askey-scheme 
of  hypergeometric  polynomials  [1].  In  this  section,  we  briefly  review  the  theory  of 
hypergeometric  orthogonal  polynomials.  We  adopt  the  notation  of  [14]  and  [17]. 

2.1.  The  Generalized  Hypergeometric  Series.  We  first  introduce  the  Pochham- 
mer  symbol  (a)„  defined  by 


(a)n  = 


1, 


if  n  =  0, 


a(a  -I-  1)  •  •  •  (a  -I-  n  —  1),  if  n  =  1, 2, 3, . . . . 
In  terms  of  Gamma  function,  we  have 

r(a  -I-  n) 


(ll)n  — 


n  >  0. 


r(a)  ’ 

The  generalized  hypergeometric  series  rFg  is  defined  by 

17  /„  „  .r  u  .  (ai)fe  •  •  •  (Or)fc 

rTs(ai)  •  ■  ■  J  O-r!  bi,  -  ■  ■  ,bs',  Z)  —  y  rrp-  — , 

^  (bi)k---(bs)k  kl 


(2.1) 


(2.2) 


(2.3) 


where  bi  yf  0,  —1,  —2, . . .  for  i  =  {1, . . . ,  s}  to  ensure  the  denominator  factors  in  the 
terms  of  the  series  are  never  zero.  Clearly,  the  ordering  of  the  numerator  parameters 
and  of  the  denominator  parameters  are  immaterial.  The  radius  of  convergence  p  of 
the  hypergeometric  series  is 

{oo  if  r  <  s  -I-  1, 

1  if  r  =  s  -I-  1,  (2.4) 

0  if  r  >  s  -I-  1. 


Some  elementary  cases  of  the  hypergeometric  series  are: 

•  Exponential  series  oTb, 

•  Binomial  series  iFq, 

•  Gauss  hypergeometric  series  27^1- 

If  one  of  the  numerator  parameters  ai,i  =  l,...,r  is  a  negative  integer,  say 
oi  =  — n,  the  hypergeometric  series  (2.3)  terminates  at  the  n*^-term  and  becomes  a 
polynomial  in  z, 


rFg  (  11-5  ‘  5  Hr  5  ‘  ‘  ‘  ^  bgj  z')  —  y  ] 

k=0 


(-n)fc  ■  ■  ■  {ar)k  z'^ 

{bi)k- ■  ■  {bs)k  kl' 


(2.5) 
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2.2.  Properties  of  the  Orthogonal  Polynomials.  A  system  of  polynomials 
{Qn{x),n  G  Af}  where  Qn(x)  is  a  polynomial  of  exact  degree  n  and  Af  =  {0, 1,  2, . . .  } 
or  Af  =  {0,1,...,  A^}  for  a  finite  nonnegative  integer  N,  is  an  orthogonal  system  of 
polynomials  with  respect  to  some  real  positive  measure  (p  if  the  following  orthogonality 
relations  are  satisfied 

/  Qn{x)Qm{x)d(l){x)  =  hlSnm,  (2.6) 

Js 

where  S  is  the  support  of  the  measure  (p  and  the  are  nonzero  constants.  The 
system  is  called  orthonormal  if  =  1. 

The  measure  p  often  has  a  density  w{x)  or  weights  w{i)  at  points  Xi  in  the  discrete 
case.  The  relations  (2.6)  then  become 


Qn{x)Qm{x)w{x)dx  =  n,m  &  Af , 

JS 


in  the  continuous  case,  or 


M 

^  ^  Qn{,Xi)Qrn{,Xi)w{^Xpj  —  ^  G  Af ,  (^-S) 

i=0 

in  the  discrete  case  where  it  is  possible  that  M  =  oo. 

The  density  w{x),  or  weights  w(i)  in  the  discrete  case,  is  also  commonly  referred 
as  the  weighting  function  in  the  theory  of  orthogonal  polynomials.  It  will  be  shown 
later  that  the  weighting  functions  for  some  orthogonal  polynomials  are  identical  to 
certain  probability  functions.  For  example,  the  weighting  function  for  the  Hermite 
polynomials  is  the  same  as  probability  density  function  of  the  Gaussian  random  vari¬ 
ables.  This  fact  plays  an  important  role  in  representing  stochastic  processes  with 
orthogonal  polynomials. 

All  orthogonal  polynomials  {Qn{x)}  satisfy  a  three-term  recurrence  relation 

xQn{x^  —  A^(5yj^l(x)  (Ajj  Cri}Qn(^X^  CnQn—l{x^<^  U  ^  1,  (^A) 

where  A„,C„  yf  0  and  >  0.  Together  with  Q-i{x)  =  0  and  Qo{x)  =  1,  all 

Qn{x)  can  be  determined  by  the  recurrence  relation. 

It  is  well  known  that  continuous  orthogonal  polynomials  satisfy  the  second-order 
differential  equation 

s{x)y”  +  T{x)y' +  \y  =  (2.10) 

where  s{x)  and  t(x)  are  polynomials  of  at  most  second  and  first  degree,  respectively, 
and 

A  =  A„  =  —nr'  —  -n(n  —  l)s"  (2.11) 

are  the  eigenvalues  of  the  differential  equation;  the  orthogonal  polynomials  y{x)  = 
yn{x)  are  the  eigenfunctions. 

In  the  discrete  case,  we  introduce  the  forward  and  backward  difference  operator, 
respectively 


d^f{x)  =  f{x  +  l)- f{x)  and  V/(x)  = /(x)  - /(x  -  1). 


(2.12) 
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The  difference  equation  corresponding  to  the  differential  equation  (2.10)  is 

s{x)AVy{x)  +  T{x)Ay{x)  +  \y{x)  =  0.  (2-13) 

Again  s(x)  and  t(x)  are  polynomials  of  at  most  second  and  first  degree,  respectively; 
X  =  Xn  are  eigenvalues  of  the  difference  equation,  and  the  orthogonal  polynomials 
y{x)  =  yn{x)  are  the  eigenfunctions. 

All  orthogonal  polynomials  can  be  obtained  by  repeatedly  applying  the  differential 
operator  as  follows 


1  d” 

Qn{x)  =  [u;(a:)s"(x)] .  (2.14) 

w[x)  dx^ 

In  the  discrete  case,  the  differential  operator  (d/dx)  is  replaced  by  the  backward  dif¬ 
ference  operator  V.  A  constant  factor  can  be  introduced  for  normalization.  Equation 
(2.14)  is  referred  as  the  generalized  Rodriguez  formula,  named  after  J.  Rodriguez  who 
first  discovered  the  specific  formula  for  Legendre  polynomials. 

2.3.  The  Askey-scheme.  The  Askey-scheme,  which  can  be  represented  as  a 
tree  structure  shown  in  figure  2.1,  classifies  the  hypergeometric  orthogonal  polynomi¬ 
als  and  indicates  the  limit  relations  between  them.  The  ‘tree’  starts  with  the  Wilson 
polynomials  and  the  Racah  polynomials  on  the  top.  They  both  belong  to  the  class 
4T3  of  the  hypergeometric  orthogonal  polynomials  (2.5).  The  Wilson  polynomials  are 
continuous  polynomials  and  the  Racah  polynomials  are  discrete.  The  lines  connecting 
different  polynomials  denote  the  limit  transition  relationships  between  them,  which 
imply  that  polynomials  at  the  lower  end  of  the  lines  can  be  obtained  by  taking  the 
limit  of  one  parameter  from  their  counterparts  on  the  upper  end.  For  example,  the 
limit  relation  between  Jacobi  polynomials  Pn°‘’^\x)  and  Hermite  polynomials  Hffx) 
is 


lim 


Hnjx) 

2”n!  ’ 


and  between  Meixner  polynomials  M„{x;  P,c)  and  Charlier  polynomials  C„(x;a)  is 

lim  (  X]  f3,  °  ^  )  =  Cn{x-,  a). 

0^00  \  a  +  p  J 

For  a  detailed  account  of  the  limit  relations  of  Askey-scheme,  the  interested  reader 
should  consult  [14]  and  [17]. 

The  orthogonal  polynomials  associated  with  the  Wiener- Askey  polynomials  chaos 
include:  Hermite,  Laguerre,  Jacobi,  Charlier,  Meixner,  Krawtchouk  and  Hahn  poly¬ 
nomials.  A  survey  with  their  definitions  and  properties  can  be  found  in  the  appendix 
of  this  paper. 

3.  The  Original  Wiener  Polynomial  Chaos.  The  Homogeneous  Chaos  ex¬ 
pansion  was  first  proposed  by  Wiener  [19];  it  employs  the  Hermite  polynomials  in 
terms  of  Gaussian  random  variables.  According  to  the  theorem  by  Cameron  &  Mar¬ 
tin  [3],  it  can  approximate  any  functionals  in  L2{C)  and  converges  in  the  L2(C) 
sense.  Therefore,  Hermite-Chaos  provides  a  means  for  expanding  second-order  ran¬ 
dom  processes  in  terms  of  orthogonal  polynomials.  Second-order  random  processes 
are  processes  with  finite  variance,  and  this  applies  to  most  physical  processes.  Thus, 
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Meixner 

2^1(2)  -  Jacobi  Meixner  Krawtchouk 

Pollaczek 


Fig.  2.1.  The  Askey  scheme  of  orthogonal  polynomials 


a  general  second-order  random  process  X(9),  viewed  as  a  function  of  9  as  the  random 
event,  can  be  represented  in  the  form 

X{9)  =  aoHo 

00 

+  ^  ai^Hi{^i^{9)) 

00  ii 

ii—1  i2  —  l 

00  ii  i2 

+  EEE  ^41^2^3  -^3  (en(0),62(^^).63W) 

il  =  l  12  =  1  13  =  1 

+  • • •  )  (3-1) 

where  . . . ,  denotes  the  Hermite-Chaos  of  order  n  in  the  variables  {^i^, . . . , 

where  the  Hn  are  Hermite  polynomials  in  terms  of  the  standard  Gaussian  variables  ^ 
with  zero  mean  and  unit  variance.  Here  ^  denotes  the  vector  consisting  of  n  indepen¬ 
dent  Gaussian  variables  The  above  equation  is  the  discrete  version  of 

the  original  Wiener  polynomial  chaos  expansion,  where  the  continuous  integrals  are 
replaced  by  summations.  The  general  expression  of  the  polynomials  is  given  by 


(3.2) 
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For  notational  convenience,  equation  (3.1)  can  be  rewritten  as 

OO 

X(0)  =  ^a,vl/^(O,  (3.3) 

j=o 

where  there  is  a  one-to-one  correspondence  between  the  functions  , . . . ,  ^j^)  and 

'l'j(^).  The  polynomial  basis  {'I'j}  of  Hermite-Chaos  forms  a  complete  orthogonal 
basis,  i.e.. 


<  >=<  >  dij,  (3.4) 

where  6ij  is  the  Kronecker  delta  and  <  •,•  >  denotes  the  ensemble  average.  This  is 
the  inner  product  in  the  Hilbert  space  determined  by  the  support  of  the  Gaussian 
variables 


<  fiOgie  >=  I  fiOgiOwm^  (3.5) 

with  weighting  function 


w{0 


\/(27r)" 


(3.6) 


What  distinguishes  the  Hermite-Chaos  expansion  from  other  possible  expansions  is 
that  the  basis  polynomials  are  Hermite  polynomials  in  terms  of  Gaussian  variables 
and  are  orthogonal  with  respect  to  the  weighting  function  W  (^)  that  has  the  form  of 
n-dimensional  independent  Gaussian  probability  density  function. 

4.  The  Wiener- Askey  Polynomial  Chaos.  The  Hermite-Chaos  expansion 
has  been  proved  to  be  effective  in  solving  stochastic  differential  equations  with  Gaus¬ 
sian  inputs  as  well  as  certain  types  of  non-Gaussian  inputs  [9],  [8],  [7],  [20];  this  can 
be  justified  by  the  Cameron-Martin  theorem  [3] .  However,  for  general  non-Gaussian 
random  inputs,  the  optimal  exponential  convergence  rate  will  not  be  realized.  In  some 
cases  the  convergence  rate  is  in  fact  severely  deteriorated. 

In  order  to  deal  with  more  general  random  inputs,  we  introduce  the  Wiener- 
Askey  polynomial  chaos  expansion  as  a  generalization  of  the  original  Wiener-Chaos 
expansion.  The  expansion  basis  is  the  complete  polynomial  basis  from  the  Askey- 
scheme  (see  section  2.3).  Similar  to  section  3,  we  represent  the  general  second-order 
random  process  X{9)  as 

X{0)  =  aolo 

OO 

+  X/  ‘^ii-^l(Cii(^)) 

OO  ii 

ii  —  1  i2  —  l 
oo  ii  12 

+  E  E  E 

il  =  l  12  =  1  13  =  1 


(4.1) 
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where  IniQi ;  ■  •  ■ )  Cin)  denotes  the  Wiener-Askey  polynomial  chaos  of  order  n  in  terms 
of  the  random  vector  ^  . . .  ,Q^).  In  the  Wiener-Askey  chaos  expansion,  the 

polynomials  In  are  not  restricted  to  Hermite  polynomials  but  rather  can  be  all  types  of 
the  orthogonal  polynomials  from  the  Askey-scheme  in  figure  2.1.  Again  for  notational 
convenience,  we  rewrite  equation  (4.1)  as 

OO 

A(0)  =  ^c,<I>,(C),  (4.2) 

3=0 

where  there  is  a  one-to-one  correspondence  between  the  functions  In{Cin  •  ■  •  )Ci„)  ^^id 
<I)j  (^).  Since  each  type  of  polynomials  from  the  Askey-scheme  form  a  complete  basis 
in  the  Hilbert  space  determined  by  their  corresponding  support,  we  can  expect  each 
type  of  Wiener-Askey  expansion  to  converge  to  any  L2  functional  in  the  L2  sense  in 
the  corresponding  Hilbert  functional  space  as  a  generalized  result  of  Cameron-Martin 
theorem  ([3]  and  [16]).  The  orthogonality  relation  of  the  Wiener-Askey  polynomial 
chaos  takes  the  form 

<  >=<  >  S,j,  (4.3) 

where  dij  is  the  Kronecker  delta  and  <  • ,  •  >  denotes  the  ensemble  average  which  is 
the  inner  product  in  the  Hilbert  space  of  the  variables  C 

<  fiOgiC)  >=  J  f{C)9{C)w{C)dC,  (4.4) 

or 

<  nOgiC)  >=  E  fiOgiOwiO  (4.5) 

c 

in  the  discrete  case.  Here  H^(C)  is  the  weighting  function  corresponding  to  the  Wiener- 
Askey  polynomials  chaos  basis  {4)^}:  see  the  appendix  for  detailed  formulas. 

As  pointed  out  in  appendix,  some  types  of  orthogonal  polynomials  from  the  Askey- 
scheme  have  weighting  functions  the  same  as  the  probability  function  of  certain  types 
of  random  distributions.  In  practice,  we  then  choose  the  type  of  independent  vari¬ 
ables  C  in  the  polynomials  {4>i(C)}  according  to  the  type  of  random  distributions  as 
shown  in  table  4.1.  It  is  clear  that  the  original  Wiener  polynomial  chaos  corresponds 


Random  variables  C 

Wiener-Askey  chaos  {4>(C)} 

Support 

Continuous 

Gaussian 

Hermite-Chaos 

(  — OO,  oo) 

Gamma 

Laguerre-Chaos 

[0,oo) 

Beta 

Jacobi-Chaos 

[a,  6] 

Uniform 

Legendre-Chaos 

[a,b] 

Discrete 

Poisson 

Charlier-Chaos 

{0,1,2,...} 

Binomial 

Krawtchouk-Chaos 

(0,1,..., iV} 

Negative  Binomial 

Meixner-Chaos 

{0,1,2,...} 

Hypergeometric 

Hahn-Chaos 

{0,1,..., iV} 

Table  4.1 

The  correspondence  of  the  type  of  Wiener-Askey  polynomial  chaos  and  their  underlying  random 
variables  (N  >  0  is  a  finite  integer). 


to  the  Hermite-Chaos  and  is  a  subset  of  the  Wiener-Askey  polynomial  chaos.  The 
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Hermite-,  Laguerre-  and  Jacobi-Chaos  are  continuous  chaos,  while  Charlier-,  Meixner- 
,  Krawtchouk-  and  Hahn-Chaos  are  discrete  chaos.  It  is  worthy  mentioning  that  the 
Legendre  polynomials,  which  is  a  special  case  of  the  Jacobi  polynomials  with  param¬ 
eters  a  =  P  =  0  (section  A. 1.3),  corresponds  to  an  important  distribution  —  the 
uniform  distribution.  Due  to  the  importance  of  the  uniform  distribution,  we  list  it 
separately  in  the  table  and  term  the  corresponding  chaos  expansion  as  the  Legendre- 
Chaos. 

5.  Applications  of  Wiener-Askey  Polynomial  Chaos.  In  this  section  we 
apply  the  Wiener-Askey  polynomial  chaos  to  solution  of  stochastic  differential  equa¬ 
tions.  We  first  introduce  the  general  procedure  of  applying  the  Wiener-Askey  poly¬ 
nomial  chaos,  then  we  solve  a  specific  stochastic  ordinary  differential  equation  with 
different  types  of  random  inputs.  We  demonstrate  the  convergence  rates  of  Wiener- 
Askey  expansion  by  comparing  the  numerical  results  with  the  corresponding  exact 
solution. 

5.1.  General  Procedure.  Let  us  consider  the  stochastic  differential  equation 

/:(x,t,6»;M)  = /(x,t;0),  (5.1) 

where  u  :=  u{y:,t;9)  is  the  solution  and  f{x,t;9)  is  the  source  term.  Operator  £ 
generally  involves  differentiations  in  space/time  and  can  be  nonlinear.  Appropriate 
initial  and  boundary  conditions  are  assumed.  The  existence  of  random  parameter  9  is 
due  to  the  introduction  of  uncertainty  into  the  system  via  boundary  conditions,  initial 
conditions,  material  properties,  etc.  The  solution  u,  which  is  regarded  as  a  random 
process,  can  be  expanded  by  the  Wiener-Askey  polynomial  chaos  as 

p 

■u(x,t;6»)  =  ^■u*(x,t)d>,(C(6')).  (5.2) 

i=0 

Note  here  the  infinite  summation  has  been  truncated  at  the  finite  term  P.  The  above 
representation  can  be  considered  as  a  spectral  expansion  in  the  random  dimension  9, 
and  the  random  trial  basis  {d)^}  is  the  Askey-scheme  based  orthogonal  polynomials 
discussed  in  section  4.  The  total  number  of  expansion  terms  is  {P  +  1),  and  is  de¬ 
termined  by  the  dimension  (n)  of  random  variable  C  ^nd  the  highest  order  (p)  of  the 
polynomials  {4)^}; 


(P+l)  =  b^  (5.3) 

nlpl 

Upon  substituting  equation  (5.2)  into  the  governing  equation  (5.1),  we  obtain 

^  t,  9;  ^  =  /(x,  t;  9).  (5.4) 

A  Galerkin  projection  of  the  above  equation  onto  each  polynomial  basis  {4)^}  is  then 
conducted  in  order  to  ensure  the  error  is  orthogonal  to  the  functional  space  spanned 
by  the  finite-dimensional  basis  {4*^}, 


£ 


t,9;  >  u^^i  ,4>fc  )  = 


i=0 


(/,4>fe),  k  =  0,l,---,P- 


(5.5) 
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By  using  the  orthogonality  of  the  polynomial  basis,  we  can  obtain  a  set  of  (P  +  1) 
coupled  equations  for  each  random  mode  Ui{x,  t)  where  z  =  {0, 1, . . . ,  P}.  It  should 
be  noted  that  by  utilizing  the  Wiener-Askey  polynomial  chaos  expansion  (5.2),  the 
randomness  is  effectively  transferred  into  the  basis  polynomials.  Thus,  the  gov¬ 
erning  equations  for  the  expansion  coefficients  Ui  resulted  from  equation  (5.5)  are 
deterministic.  Discretizations  in  space  x  and  time  t  can  be  carried  out  by  any  conven¬ 
tional  deterministic  techniques,  e.g.,  Runge-Kutta  solvers  in  time  and  the  spectral//ip 
element  method  in  space  for  high  accurate  solution  in  complex  geometry  [13]. 

5.2.  Stochastic  Ordinary  Differential  Equation.  We  consider  the  ordinary 
differential  equation 


=  -ky,  y{0)  =  y,  (5.6) 

where  the  decay  rate  coefficient  k  is  considered  to  be  a  random  variable  k{0)  with 
certain  distribution  and  mean  value  k.  The  probability  function  is  f{k)  for  the  con¬ 
tinuous  case  or  /(fcz)  for  the  discrete  case.  The  deterministic  solution  is 

y{t)  =  yoe~’"*  (5.7) 


and  the  mean  of  stochastic  solution  is 


y{t)=yj  e  ^*f{k)dk  or  y{t)=y'^e  (5.8) 


corresponding  to  the  continuous  and  discrete  distributions,  respectively.  The  inte¬ 
gration  and  summation  are  taken  within  the  support  defined  by  the  corresponding 
distribution. 

By  applying  the  Wiener-Askey  polynomial  chaos  expansion  (4.2)  to  the  solution 
y  and  random  input  k 

p  p 

y{t)  =^yi{t)^i,  k  =  ^ki^,  (5.9) 

i— 0  i— 0 

and  substituting  the  expansions  into  the  governing  equation,  we  obtain 


p 


E 


dyt{t)^ 

dt  * 


p  p 


jkiyj{t). 

i=o  i=o 


(5.10) 


We  then  project  the  above  equation  onto  the  random  space  spanned  by  the  orthogonal 
polynomial  basis  {4>i}  by  taking  the  inner  product  of  the  equation  with  each  basis. 
By  taking  <  .,<!>;  >  and  utilizing  the  orthogonality  condition  (4.3),  we  obtain  the 
following  set  of  equations: 


dyi{t) 

dt 


1 

<W> 


P  P 

EE  (^)  j 

i—0  j—0 


(5.11) 


where  =<  >.  Note  that  the  coefficients  are  smooth  and  thus  any  standard 

ODE  solver  can  be  employed  here.  In  the  following  the  standard  second-order  Runge- 
Kutta  scheme  is  used. 
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5.3.  Numerical  Results.  In  this  section  we  present  numerical  results  of  the 
stochastic  ordinary  differential  equation  by  the  Wiener-Askey  polynomial  chaos  ex¬ 
pansion.  For  the  purpose  of  benchmarking,  we  will  arbitrarily  assume  the  type  of 
distributions  of  the  decay  parameter  k  and  employ  the  corresponding  Wiener-Askey 
chaos  expansion,  although  in  practice  there  is  certainly  more  favorable  assumptions 
about  k  depending  on  the  specific  physical  background.  We  define  the  two  error 
measures  for  the  mean  and  variance  of  the  solution 


.{t)  = 


y{^)  2/exact  (t) 


2/exact  (t) 


^var(t) 


rrexact(t) 
^exact  (t) 


(5.12) 


where  y{t)  =  E[y{t)]  is  the  mean  value  of  y{t)  and  a{t) 


E 


iyit)-yit)y 


is  the 


variance  of  the  solution.  The  initial  condition  is  fixed  to  be  y  =  1  and  the  integration 
is  performed  up  to  t  =  1  (nondimensional  time  units). 


5.3.1.  Gaussian  Distribution  and  Hermite-Chaos.  In  this  section  the  dis¬ 
tribution  of  k  is  assumed  to  be  a  Gaussian  random  variable  with  probability  density 
function 


m 


(5.13) 


which  has  zero  mean  value  (fc  =  0)  and  unit  variance  (ct^  =  1).  The  exact  stochastic 
mean  solution  is 


y{t)  = 


(5.14) 


The  Hermite-Chaos  from  the  Wiener-Askey  polynomial  chaos  family  is  employed  as 
a  natural  choice  due  to  the  fact  that  the  random  input  is  Gaussian.  Figure  5.1  shows 


2 

1.5 


0.5 

c 

_o 

o 

CO 

-0.5 


■1.5 


Vo  (mean) 

y, 

72 


Vs 

Va 

Deterministic 


-2 


0 


'  o'5  ' 

Time 


Fig.  5.1.  Solution  with  Gaussian  random  input  by  4th-order  Hermite-Chaos;  Left:  Solution  of 
each  random  mode,  Right:  Error  convergence  of  the  mean  and  the  variance. 


the  solution  by  the  Hermite-Chaos  expansion.  The  convergence  of  errors  of  the  mean 
and  variance  as  the  number  of  expansion  terms  increases  is  shown  on  semi-log  plot, 
and  it  is  seen  that  the  exponential  convergence  rate  is  achieved.  It  is  also  noticed  that 
the  deterministic  solution  remains  constant  as  the  mean  value  of  k  is  zero;  however 
the  mean  of  the  stochastic  solution  (random  mode  with  index  0,  yo)  is  nonzero  and 
grows  with  time. 
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5.3.2.  Gamma  Distribution  and  Laguerre- Chaos.  In  this  section  we  as¬ 
sume  the  distribution  of  the  decay  parameter  k  is  the  Gamma  distribution  with  PDF 
of  the  form 


0<fc<oo,  a>-l.  (5.15) 

r(Q;  -I-  1) 

The  mean  and  variance  of  k  are:  fik  =  k  =  a  +  1  and  cr^  =  a  -I-  1,  respectively.  The 
mean  of  stochastic  solution  is 


yW  =^(l+^^)a+i-  (5-16) 

The  special  case  of  a  =  0  corresponds  to  another  important  distribution:  the  exponential 
distribution.  Because  the  random  input  has  a  Gamma  distribution,  we  employ  the 


Fig.  5.2.  Solution  with  Gamma  random  input  by  4th-order  Laguerre- Chaos;  Left:  Solution  of 
each  mode  (a  =  0;  exponential  distribution),  Right:  Error  convergence  of  the  mean  and  the  variance 
with  different  a. 


Laguerre-Ghaos  as  the  specific  Wiener- Askey  chaos  (see  table  4.1).  Figure  5.2  shows 
the  evolution  of  each  solution  mode  over  time,  together  with  the  convergence  of  the 
errors  of  the  mean  and  the  variance  with  different  values  of  parameter  a.  The  special 
case  of  exponential  distribution  is  included  {a  =  0).  Again  the  mean  of  stochastic 
solution  and  deterministic  solution  show  significant  difference.  As  a  becomes  larger, 
the  spread  of  the  Gamma  distribution  is  larger  and  this  leads  to  larger  errors  with 
fixed  number  of  Laguerre-Ghaos  expansion.  However,  the  exponential  convergence 
rate  is  still  realized. 

5.3.3.  Beta  Distribution  and  Jacobi-Chaos.  We  now  assume  the  distribu¬ 
tion  of  the  random  variable  k  to  be  the  Beta  distribution  with  probability  density 
function  of  the  form 


f{k]a,P) 


{l-k)°^{l  +  kf 
2“+/5+iH(a-Fl,/3-Fl)’ 


—  1  <  fc  <  1,  a,P  >  — 1, 


(5.17) 


where  B{a,j3)  is  the  Beta  function  defined  as  B{p,q)  =  F(p)F((7)/F(p -|-  q).  We  then 
employ  the  Jacobi-Ghaos  expansion  which  has  the  weighting  function  in  the  form  of 
the  Beta  distribution.  An  important  special  case  is  a  =  /3  =  0  when  the  distribution 
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Fig.  5.3.  Solution  with  Beta  random  input  by  4th-order  Jacobi-Chaos;  Left:  Solution  of  each 
mode  (a  =  f}  =  0;  Legendre- Chaos),  Right:  Error  convergence  of  the  mean  and  the  variance  with 
different  a  and  fS. 


becomes  the  uniform  distribution  and  the  corresponding  Jacobi-Chaos  becomes  the 
Legendre-Chaos. 

Figure  5.3  shows  the  solution  by  the  Jacobi-Chaos.  On  the  left  is  the  evolution 
of  all  random  modes  of  the  Legendre-Chaos  {a  =  f3  =  Q)  with  uniformly  distributed 
random  input.  In  this  case,  k  has  zero  mean  value  and  the  deterministic  solution 
remains  constant,  but  the  mean  of  stochastic  solution  grows  over  time.  The  con¬ 
vergence  of  errors  of  the  mean  and  the  variance  of  the  solution  with  respect  to  the 
order  of  Jacobi-Chaos  expansion  is  shown  on  the  semi- log  scale,  and  the  exponential 
convergence  rate  is  obtained  with  different  sets  of  parameter  values  a  and  /J. 

5.3.4.  Poisson  Distribution  and  Char  Her- Chaos.  We  now  assume  the  dis¬ 
tribution  of  the  decay  parameter  k  to  be  Poisson  of  the  form 

\  ^ 

/(fc;  A)  =  e-^  — ,  fc  =  0, 1, 2, . . . ,  A  >  0.  (5.18) 

The  mean  and  variance  of  k  are:  Hk  =  k  =  X  and  =  A,  respectively.  The  analytic 
solution  of  the  mean  stochastic  solution  is 

y{t)  =  (5.19) 

The  Charlier-Chaos  expansion  is  employed  to  represent  the  solution  process  and  the 
results  with  fourth-order  expansion  are  shown  in  figure  5.4.  Once  again  we  see  the 
noticeable  difference  between  the  deterministic  solution  and  the  mean  of  stochastic 
solution.  Exponential  convergence  rate  is  obtained  for  different  values  of  parameter 

A. 

5.3.5.  Binomial  Distribution  and  Krawtchouk-Chaos.  In  this  section  the 
distribution  of  the  random  input  k  is  assumed  to  be  binomial 

/(A:;p,iV)=  (^^^/(l-p)'^-^  0<p<  1,  A:  =  0,l,...,iV.  (5.20) 

The  exact  mean  solution  of  (5.6)  is 

y{t)  =y[l-  (l-e"‘)p]^. 


(5.21) 
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Fig.  5.4.  Solution  with  Poisson  random  input  by  4th-order  Charlier- Chaos;  Left:  Solution  of 
each  mode  (\  =  1),  Right:  Error  convergence  of  the  mean  and  the  variance  with  different  A. 


Fig.  5.5.  Solution  with  binomial  random  input  by  fth-order  Krawtchouk-Chaos;  Left:  Solution 
of  each  mode  (p  =  0.5,  A"  =  b)),  Right:  Error  convergence  of  the  mean  and  the  variance  with 
different  p  and  N. 


Figure  5.5  shows  the  solution  with  4th-order  Krawtchouk-Chaos.  With  different 
parameter  sets,  Krawtchouk-Chaos  expansion  correctly  approximates  the  exact  solu¬ 
tion,  and  the  convergence  rate  with  respect  to  the  order  of  expansion  is  exponential. 

5.3.6.  Negative  Binomial  Distribution  and  Meixner- Chaos.  In  this  sec¬ 
tion  we  assume  the  distribution  of  the  random  input  of  k  is  the  negative  binomial 
distribution 

/(A:;/3,c)  =  ^(l-c)^c^  0<c<  1,  /3>0,  A:  =  0,1,....  (5.22) 

k\ 

In  case  of  (3  being  integer,  it  is  often  called  the  Pascal  distribution.  The  exact  mean 
solution  of  (5.6)  is 

y{t)  =  y  I  ^  ■  (5-23) 

The  Meixner-Chaos  is  chosen  since  the  random  input  is  negative  binomial  (see 
table  4.1).  Figure  5.6  shows  the  solution  with  4th-order  Meixner-Chaos.  Exponential 
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Yg  (mean) 


Fig.  5.6.  Solution  with  negative  binomial  random  input  by  Jfth-order  Meixner-Chaos;  Left: 
Solution  of  each  mode  (0  =  l,c  =  0.5)),  Right:  Error  convergence  of  the  mean  and  the  variance 
with  different  0  and  c. 


convergence  rate  is  observed  by  the  Meixner-Chaos  approximation  with  different  sets 
of  parameter  values. 

5.3.7.  Hypergeometric  Distribution  and  Hahn-Chaos.  We  now  assume 
the  distribution  of  the  random  input  k  is  hypergeometric 


f{k;a,P,N)=  >  fc  =  0,l,...,iV,  a,/3>iV.  (5.24) 

\  N  ) 


Fig.  5.7.  Solution  with  hypergeometric  random  input  by  4-lh-order  Hahn-Chaos;  Left:  Solution 
of  each  mode  (a  =  /3  =  5,  =  4.)),  Right:  Error  convergence  of  the  mean  and  the  variance  with 

different  a,  /3  and  N. 


In  this  case,  the  optimal  Wiener- Askey  polynomial  chaos  is  the  Hahn-Chaos  (table 
4.1).  Figure  5.7  shows  the  solution  by  4th-order  Hahn-Chaos.  It  can  be  seen  from  the 
semi-log  plot  of  the  errors  of  the  mean  and  variance  of  the  solution  that  exponential 
convergence  rate  is  obtained  with  respect  to  the  order  of  Hahn-Chaos  expansion  for 
different  sets  of  parameter  values. 
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5.4.  Efficiency  of  Wiener- Askey  Chaos  Expansion.  We  have  demonstrated 
the  exponential  convergence  of  the  Wiener-Askey  polynomial  chaos  expansion.  From 
the  results  above,  we  notice  that  it  normally  takes  an  expansion  order  P  =  2  ~  4  for 
the  error  of  the  mean  solution  to  reach  the  order  of  0(10“^).  Equation  (5.11)  shows 
that  the  Wiener-Askey  chaos  expansion  with  highest  order  of  P  results  in  a  set  of 
{P  +  1)  coupled  ODEs.  Thus,  the  computational  cost  is  slightly  more  than  {P  +  1) 
times  of  that  of  a  single  realization  of  the  deterministic  integration.  On  the  other 
hand,  if  the  Monte-Carlo  simulation  is  used,  it  normally  requires  0(10"^)  ~  0(10®) 
number  of  realizations  to  reduce  the  error  of  the  mean  solution  to  0(10“®).  For 
example,  if  k  is  an  exponentially  distributed  random  variable,  the  error  convergence 
of  the  mean  solution  of  the  Monte-Carlo  simulation  is  shown  in  table  5.1. 


N 

1  X  10^ 

1  X  10® 

1  X  lO"^ 

1  X  10® 

^mean 

4.0  X  10“^ 

1.1  X  10“^ 

5.1  X  10“® 

6.5  X  10“4 

Table  5.1 


Error  convergence  of  the  mean  solution  by  Monte-Carlo  simulation:  N  is  the  number  of  real¬ 
izations  and  Emean  is  the  error  of  mean  solution  defined  in  (5,12);  Random  input  has  exponential 
distribution. 


Monte-Carlo  simulations  with  other  types  of  random  inputs  as  discussed  in  this 
paper  have  also  been  conducted  and  the  results  are  similar.  The  actual  numerical 
values  of  the  errors  with  given  number  of  realizations  may  vary  depending  on  the 
property  of  random  number  generators  used,  but  the  order  of  magnitude  should  be 
the  same.  Techniques  such  as  variance  reduction  are  not  used.  Although  such  tech¬ 
niques,  if  applicable,  can  greatly  speed  up  Monte-Carlo  simulation  by  an  order  or  more 
depending  on  the  specific  problem,  the  advantage  of  Wiener-Askey  polynomial  chaos 
expansion  is  obvious.  For  the  ordinary  differential  equation  discussed  in  this  paper, 
speed-up  of  order  0(10®)  ~  O(IO^)  compared  with  straight  Monte-Carlo  simulations 
can  be  expected.  However,  for  more  complicated  problems  where  there  exist  multi¬ 
dimensional  random  inputs,  the  multi-dimensional  Wiener-Askey  chaos  is  needed. 
The  total  number  of  expansion  terms  increases  fast  for  large  dimensional  problems 
(see  equation(5.3)).  Thus  the  efficiency  of  the  chaos  expansion  can  be  reduced. 

6.  Representation  of  Arbitrary  Random  Inputs.  As  demonstrated  above, 
with  appropriately  chosen  Wiener-Askey  polynomial  chaos  expansion  according  to  the 
type  of  the  random  input,  optimal  exponential  convergence  rate  of  the  chaos  expansion 
can  be  realized.  In  practice,  we  often  encounter  distributions  of  random  inputs  not 
belonging  to  the  basic  types  of  distributions  listed  in  table  4.1,  or  even  when  they 
do  belong  to  certain  basic  types,  the  correspondence  is  not  explicitly  known.  In  this 
case,  we  need  to  project  the  input  process  onto  the  Wiener-Askey  polynomial  chaos 
basis  directly  in  order  to  solve  the  differential  equation. 

Let  us  assume  in  the  stochastic  ODE  of  equation  (5.6)  that  the  distribution  of 
the  decay  parameter  k  is  known  in  the  form  of  probability  function  /(fc).  The  repre¬ 
sentation  of  k  by  the  Wiener-Askey  polynomial  chaos  expansion  takes  the  form 


p 

k  =  '^  ki^i, 

i=0 


ki  = 


<  > 

<  >  ’ 


(6.1) 


where  the  operation  <  •,  •  >  denotes  the  inner  product  in  the  Hilbert  space  spanned 
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by  the  Wiener-Askey  chaos  basis  i.e., 

/'fe^*(C)5(C)c?C.  or  k,=  (6.2) 

I  I  j 

where  g{x)  and  g{xi)  is  the  probability  function  of  the  random  variable  C  in  the 
Wiener-Askey  polynomial  chaos  for  continuous  and  discrete  cases,  respectively.  The 
underlying  assumption  here  is  that  the  random  variable  C  is  fully  dependent  on  the 
target  random  variable  k.  We  notice  that  the  above  equations  are  mathematically 
meaningless  due  to  the  fact  that  the  support  of  k  and  C  are  likely  to  be  different.  In 
other  words,  the  random  variables  k  and  C  could  belong  to  two  different  probability 
spaces  {n,A,P)  with  different  event  space  ft,  cr-algebra  A  and  probability  measure 
P. 

6.1.  Analytical  Approach.  In  order  to  conduct  the  above  projection,  we  need 
to  transform  the  fully  correlated  random  variables  k  and  C  to  the  same  probability 
space.  Under  the  theory  of  probability  this  is  always  possible.  In  practice,  it  is 
convenient  to  transform  them  to  the  uniformly  distributed  probability  space  u  C 
C/(0, 1).  In  fact,  the  inverse  procedure  is  an  important  technique  for  random  number 
generation  where  one  first  generates  the  uniformly  distributed  numbers  as  the  seeds 
and  then  performs  the  inverse  transformation  according  to  the  desired  distribution 
function.  Without  loss  of  generality,  we  discuss  in  detail  the  case  when  k  and  C  are 
continuous  random  variables. 

Let  us  assume  that  the  random  variable  u  is  uniformly  distributed  in  (0,1)  and 
the  probability  density  functions  for  k  and  C  are  f{k)  and  (/(C))  respectively.  A 
transformation  of  variables  in  probability  space  shows  that 

du  =  f{k)dk  =  dF{k),  du  =  g{0dC  =  dG{0,  (6.3) 

where  F  and  G  are  the  distribution  function  of  k  and  C)  respectively. 


If  we  require  the  random  variables  k  and  C  to  be  transformed  to  the  same  uniformly 
distributed  random  variable  m,  we  obtain 

u  =  F{k)  =  G{0.  (6.5) 

After  inverting  the  above  equations,  we  obtain 

k  =  F~^{u)  =  h{u),  (  =  G~^{u)  =  l{u).  (6.6) 

Now  that  we  have  effectively  transformed  the  two  different  random  variables  k  and 
C  to  the  same  probability  space  defined  by  m  G  C/(0,1),  the  projection  (6.2)  can  be 
performed,  i.e., 

k^  =  J  kMOgiOdC 

1 

h{u)^i{l{u))du.  (6-7) 
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In  general  the  above  integral  can  not  be  integrated  analytically.  However,  it  can 
be  efficiently  evaluated  with  the  Gauss  quadrature  in  the  closed  domain  [0, 1]  with 
sufficient  accuracy.  The  analytical  forms  of  the  inversion  relations  (6.6)  are  known  for 
some  basic  distributions:  Gaussian,  exponential.  Beta,  etc  [6]. 

The  above  procedure  works  equally  well  for  the  discrete  distributions,  where  the 
inversion  procedure  is  slightly  modified  and  the  integral  in  (6.7)  is  replaced  by  sum¬ 
mation. 

6.2.  Numerical  Approach.  The  procedure  described  above  requires  the  dis¬ 
tribution  functions  F{k)  and  0(0  be  known  and  the  inverse  functions  F~^  and  G~^ 
exist  and  be  known  as  well.  In  practice,  these  conditions  are  not  always  satisfied. 
Often  we  only  know  the  probability  function  /(fc)  for  a  specific  problem.  The  proba¬ 
bility  function  g{(^)  is  known  from  the  choice  of  Wiener-Askey  polynomial  chaos  but 
the  inversion  is  not  always  known  either.  In  this  case,  we  can  perform  the  projec¬ 
tion  (6.2)  directly  by  the  Monte-Garlo  integration,  where  a  large  ensemble  of  random 
numbers  k  and  C  are  generated.  The  requirement  of  k  and  C  being  transformed  to  the 
same  probability  space  u  G  C/(0, 1)  by  equation  (6.5)  implies  that  each  pair  of  k  and 
C  has  to  be  generated  from  the  same  seed  of  uniformly  generated  random  number 
u  G  C/(0,1). 

6.3.  Results.  In  this  section  we  present  numerical  examples  of  representing  an 
arbitrarily  given  random  distribution.  More  specificly,  we  present  results  of  using 
Hermite-Ghaos  expansion  for  some  non-Gaussian  random  variables.  Although  in 
theory,  Hermite-Ghaos  converges  and  it  has  been  successfully  applied  to  some  non- 
Gaussian  processes  ([8],  [20]),  we  demonstrate  numerically  that  optimal  exponential 
convergence  rate  is  not  realized. 

6.3.1.  Approximation  of  Gamma  distribution  by  Hermite-Ghaos.  Let 

us  assume  that  the  decay  parameter  k  in  the  ODE  (5.6)  is  a  random  variable  with 
Gamma  distribution  (5.15).  We  consider  the  specific  case  of  a  =  0.  In  this  case  k  is 
a  random  variable  with  exponential  distribution  with  PDF  of  the  form 

/(A:)  =  e-^  fc>0.  (6.8) 

The  inverse  of  its  distribution  function  F{k)  (equation  (6.6))  is  known  as 

h{u)  =  F~^  (u)  = —ln{l  —  u),  uGU{0,1).  (6.9) 


We  then  use  Hermite-Ghaos  to  represent  k  instead  of  the  optimal  Laguerre-Ghaos. 
The  random  variable  C  in  equation  (6.7)  is  a  standard  Gaussian  variable  with  PDF 
g(C)  =  The  inverse  of  the  Gaussian  distribution  G(C)  is  known  as 


l{u) 


G  ^(u)  =  sign 


Co  +  Cit  +  C2t^  \ 
l  +  dit  +  d2F  +  dsF  )  ’ 


(6.10) 


where 


t  = 


In  [min(u,  1 


u)] 


2 


and 


Co  =  2.515517,  Cl  =  0.802853,  C2  =  0.010328, 
di  =  1.432788,  da  =  0.189269,  ds  =  0.001308. 
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Fig.  6.1.  Approximation  of  exponential  distribution  with  Hermite- Chaos;  Left:  The  expansion 
coefficients,  Right:  The  PDF  of  different  orders  of  approximations. 


The  formula  is  from  Hastings  [10]  and  the  numeric  values  of  the  constants  have  abso¬ 
lute  error  less  than  4.5  x  10  ^  (also  see  [6]). 

In  figure  6.1  we  show  the  result  of  the  approximation  of  the  exponential  distri¬ 
bution  by  the  Hermite-Chaos.  The  expansion  coefficients  ki  are  shown  on  the  left, 
and  we  see  the  major  contributions  of  the  Hermite-Chaos  approximation  are  from  the 
first  three  terms.  The  PDF  of  different  orders  of  the  approximations  are  shown  on  the 
right,  together  with  the  exact  PDF  of  the  exponential  distribution.  We  notice  that 
the  third-order  approximation  gives  fairly  good  result  and  fifth-order  Hermite-Chaos 
is  very  close  to  the  exact  distribution.  The  Hermite-Chaos  does  not  approximate  the 
PDF  well  at  a;  ~  0  where  the  PDF  reaches  its  peak  at  1.  In  order  to  capture  this 
rather  sharp  region,  more  Hermite-Chaos  terms  are  needed. 

The  above  result  is  the  representation  of  the  random  input  k  for  the  ODE  of 
equation  (5.6).  If  the  optimal  Wiener- Askey  chaos  is  chosen,  in  this  case  the  Laguerre- 
Chaos,  only  one  term  is  needed  to  represent  k  exactly.  We  can  expect  if  the  Hermite- 
Chaos  is  used  to  solve  the  differential  equation  in  this  case,  the  solution  would  not 
retain  the  exponential  convergence  as  realized  by  the  Laguerre-Chaos. 

In  figure  6.2  the  errors  of  mean  solution  defined  by  equation  (5.12)  with  Laguerre- 
Chaos  and  Hermite-Chaos  to  the  ODE  of  equation  (5.6)  are  shown.  The  random 
input  of  k  has  exponential  distribution  which  implies  that  the  Laguerre-Chaos  is 
the  optimal  Wiener- Askey  polynomial  chaos.  It  is  seen  from  the  result  that  the 
exponential  convergence  rate  is  not  obtained  by  the  Hermite-Chaos  as  opposed  to  the 
Laguerre-Chaos. 

6.3.2.  Approximation  of  Beta  distribution  by  Hermite-Chaos.  We  now 

assume  the  distribution  of  k  is  Beta  distribution;  see  equation  (5.17).  We  return  to 
the  more  conventional  definition  of  Beta  distribution  in  the  domain  [0, 1] 

0<k<l.  (6.11) 

B(a  -I-  1,  p  -I-  1) 

Figure  6.3  shows  the  PDF  of  first-,  third-  and  fifth-order  Hermite-Chaos  approxi¬ 
mations  to  the  Beta  random  variable.  The  special  case  of  a  =  /3  =  0  is  the  important 
uniform  distribution.  It  can  be  seen  that  the  Hermite-Chaos  approximation  converges 
to  the  exact  solution  as  the  number  of  expansion  terms  increases.  Oscillations  are  ob- 


20  D.  XIU  AND  G.  E.  KARNIADAKIS 


P 


Fig.  6.2.  Error  convergence  of  the  mean  solution  of  the  Laguerre- Chaos  and  Hermite-Chaos  to 
stochastic  ODE  with  random  input  of  the  exponential  distribution 


A 

-  3rd-order 
-  Sth-ordsr 

_ .  .  ■  \ 

1 

1 

Fig.  6.3.  PDF  of  approximations  of  Beta  distributions  by  Hermite-Chaos;  Left:  a  =  f)  =  0, 
the  uniform  distribution,  Right:  a  =  2,  f)  =  0. 


served  near  the  corners  of  the  square.  This  is  in  analogy  with  the  Gibb’s  phenomena 
which  occurs  when  Fourier  expansions  are  used  to  approximate  functions  with  sharp 
corners.  Since  all  Wiener-Askey  polynomial  chaos  expansions  can  be  considered  as 
spectral  expansions  in  the  random  dimension,  the  oscillations  here  can  be  regarded 
as  the  stochastic  Gibb ’s  phenomena.  For  uniform  distribution,  Hermite-Chaos  does 
not  work  very  well  due  to  the  stochastic  Gibb’s  phenomena  even  when  more  higher- 
order  terms  are  added.  On  the  other  hand,  the  first-order  Jacobi-Chaos  expansion  is 
already  exact.  In  addition  to  the  exponential  convergence,  the  proper  Wiener-Askey 
basis  leads  to  dramatic  lowering  of  dimensionality  of  the  problem. 

7.  Conclusion.  We  have  proposed  a  Wiener-Askey  polynomial  chaos  expansion 
to  represent  stochastic  processes  and  further  model  the  uncertainty  in  practical  appli¬ 
cations.  The  Wiener-Askey  polynomial  chaos  can  be  regarded  as  the  generalization 
of  the  homogeneous  chaos  first  proposed  by  Wiener  in  1938.  The  original  Wiener  ex¬ 
pansion  employs  the  Hermite  polynomials  in  terms  of  Gaussian  random  variables.  In 
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the  Wiener-Askey  chaos  expansion,  the  basis  polynomials  are  those  from  the  Askey- 
scheme  of  hypergeometric  orthogonal  polynomials,  and  the  underlying  variables  are 
random  variables  chosen  according  to  the  weighting  function  of  the  polynomials.  A 
general  guideline  of  choosing  the  optimal  Wiener-Askey  polynomial  chaos  according 
to  the  random  inputs  is  given.  By  solving  a  stochastic  ordinary  differential  equa¬ 
tion,  we  demonstrate  numerically  that  the  Wiener-Askey  polynomial  chaos  exhibits 
exponential  convergence  rate.  For  any  given  type  of  random  input,  the  Wiener-Askey 
polynomial  chaos  converges  in  general,  although  the  exponential  rate  is  not  retained 
if  the  optimal  chaos  is  not  chosen.  The  Wiener-Askey  polynomial  chaos  proposed  in 
the  present  paper  can  deal  with  general  random  inputs  more  effectively  than  the  orig¬ 
inal  Wiener-Hermite  chaos.  It  can  be  extended  to  more  complex  stochastic  systems 
governed  by  partial  differential  equations  without  any  fundamental  difficulties. 

Appendix  A.  Some  Important  Orthogonal  Polynomials  in  Askey-scheme. 

In  this  section  we  briefly  review  the  definitions  and  properties  of  some  important 
orthogonal  polynomials  from  Askey  scheme,  which  are  discussed  in  this  paper  for  the 
Wiener-Askey  polynomial  chaos. 

A.l.  Continuous  Polynomials. 

A. 1.1.  Hermite  Polynomial  and  Gaussian  Distribution. 

Definition: 

H^{x)  =  {2xr  2F0  •  (A.l) 

Orthogonality: 

1  2 

Hm{x)H„{x)dx  =  2^n\Smn-  (A. 2) 

V  ^  J  —  00 

Recurrence  relation: 

Hn+i{x)  -  2xHn{x)  +  2nHr,-i{x)  =  0.  (A. 3) 

Rodriguez  formula: 

e-^^H^{x)  =  (-1)”^  (e-^^)  .  (A.4) 

The  weighting  function  is  w{x)  =  e~^  from  the  orthogonality  condition  (A. 2). 
After  rescaling  x  by  the  weighting  function  is  the  same  as  the  probability  density 
function  of  a  standard  Gaussian  random  variable  with  zero  mean  and  unit  variance. 

A.  1.2.  Laguerre  Polynomial  (x)  and  Gamma  Distribution. 
Definition: 


{a  +  1)„ 
n! 


iFi(— n;  0-1-1;  x). 


(A.5) 


Orthogonality: 


^L^^\x)L[^\x)dx  = 


r(r 


1) 


a  >  — 1. 


(A.6) 


0 
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Recurrence  relation: 

(n+  -  (2n  +  a  +  1  -  x)lI^\x)  +  (n  +  a)L^“\(x)  =  0.  (A.7) 

Rodriguez  formula: 

1 

(e-V+“).  (A.8) 

Recall  that  the  Gamma  distribution  has  the  probability  density  function 

+  (A.9) 

Despite  of  the  scale  parameter  (3  and  a  constant  factor  r(a+  1),  it  is  the  same  as  the 
weighting  function  of  Laguerre  polynomial. 

A. 1.3.  Jacobi  Polynomial  P^’^\x)  and  Beta  Distribution. 

Definition: 

P^'^\x)  =  +  ^)»  Q,  ^  1-  Q,  _|_  Ij  ^  .  (A. 10) 

n!  \  I  ) 

Orthogonality: 

J  {1  -  x)°'{l  +  x)^p^'^\x)p^°''^\x)dx  =  hlSmn,  a  > -1,  (3  > -1,  (A. 11) 


where 


2a+/3+i  r(n  +  a  +  l)r(n  +  /3+l) 

2n  +  a  +  f3+l  T{n  +  a  +  (3  +  l)n\ 


Recurrence  relation: 


=  2(n+l)(n  +  a  +  /3+l)  p(a,/3)/  8 

”  ^  ’  (2n  +  a  +  /3+l)(2n  +  a  +  /3  +  2)  ^  ’ 


_ _ 

(2?i  “h  o  -|-  /3)(2?t-  ct  jS  ‘2'j 

2{n  +  a){n  +  /?) 

(^2ti  “h  q;  -|-  /3)(2?T-  -h  O  -h  /3  -h  1) 


py-«(a.') 


Rodriguez  formula: 

(1  _  ^)«(1  +  :,)/5p(-A)(^)  =  [(1  -  x)"+“(l  +  x)"+^] 

The  Beta  distribution  has  the  probability  density  function 

...  (x  —  a)^(b  —  x)“  , 

(i,-a)°+Wi;(„+i,ff+i)-  “<=^<1’. 

where  B{p,q)  is  the  Beta  function  defined  as 

_  r(p)r(g) 


B{p,  q)  = 


T{p  +  q) 


(A.12) 

(A.13) 

(A.14) 

(A.15) 


It  is  clear  that  despite  of  a  constant  factor  the  weighting  function  of  Jacobi  polynomial 
w{x)  =  (1  — x)“(l  +  x)^  from  (A. 11)  is  the  same  as  the  probability  density  function  of 
Beta  distribution  defined  in  domain  [—1, 1].  When  a  =  P  =  0,  the  Jacobi  polynomials 
become  the  Legendre  polynomials  and  the  weighting  function  is  a  constant  which 
corresponds  to  the  important  uniform  distribution. 
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A. 2.  Discrete  Polynomials. 

A. 2.1.  Charlier  Polynomial  Cn{x;a)  and  Poisson  Distribution. 
Definition: 


Cn{x-,a)=  2F0  (^-n,-x]  5“^^'  (A. 16) 

Orthogonality: 

00 

‘^Cm{x;  a)Cn{x;  a)  =  a“”e“n!(5„i„,  a  >  0.  (A. 17) 

x—0 


Recurrence  relation: 

—xCn{x]  a)  =  aCn+i(x;  a)  —  (n  +  a)C„{x;  a)  +  nCn-i{x;  a).  (A. 18) 

Rodriguez  formula: 

=  (A,19) 

where  V  is  the  backward  difference  operator  (2.12). 

The  probability  function  of  Poisson  distribution  is 

/(a:;a)  =  e-“J,  fc  =  0,l,2,....  (A.20) 

Despite  of  a  constant  factor  e““,  it  is  the  same  as  the  weighting  function  of  Charlier 
polynomials. 

A. 2. 2.  Krawtchouk  Polynomial  Kn{x;p,  N)  and  Binomial  Distribution. 
Definition: 


Kn{x;  p,  N)  =  2F1  -n,  -x;  -N;  -  ,  n  =  0, 1, . . . ,  A. 

p; 


(A.21) 


Ort  hogonality : 

N 


^  '"Km{x;p,N)Kn{x;p,N)  =  (~^)  0  <  P  <  1- 

(A.22) 


x=0 

Recurrence  relation: 


-xK{x;p,N)  =  p{N  -  n)Kn+i{x;p,  N)  -  [p{N  -  n)  +  n(l  -  p)]Kn{x;p,  N) 

+  n{l-  p)Kn-i{x;p,N).  (A. 23) 


Rodriguez  formula: 

^N\  f  p 


X  J  \l-p 


Kr,{x-,p,N)  =  V’^- 


N  —  n\  f  p 

X 


l-p 


(A.24) 


Clearly,  the  weighting  function  from  (A.22)  is  the  probability  function  of  the 
binomial  distribution. 
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A. 2. 3.  Meixner  Polynomial  c)  and  Negative  Binomial  Distri¬ 

bution. 

Definition: 


Mn{x-,  (3,  c)  =  2F1  (  -n,  -X]  /3;  1  -  - 

c 


(A.25) 


Orthogonality: 

(3,  c)Mn{x;  (3,  c)  =  ,  /3  >  0,  0  <  c  <  1.  (A. 26) 

^_n  IPjnlt  C)^^ 


-n  f 


x—0 

Recurrence  relation: 


(c  -  l)xM„(x;  (3,  c)  =  c(n  +  /3)M„+i(a;;  (3,  c)  -  [n  +  (n  +  /3)c]M„(a;;  (3,  c) 

+  nMn-i{x-,(3,c).  (A. 27) 


Rodriguez  formula: 


(/3)^ 

x! 


M„(x;/3,c)  =  V" 


(/?  +  n)a;C^ 

x! 


(A.28) 


The  weighting  function  is 

/(x)  =  ^(l-c)^c^,  0<p<l,  /3>0,  a;  =  0,1,2,....  (A.29) 

x\ 

It  can  verified  that  it  is  the  probability  function  of  negative  binomial  distribution.  In 
the  case  of  [3  being  integer,  it  is  often  called  the  Pascal  distribution. 

A. 2. 4.  Hahn  Polynomial  Qn{x;a,P,N)  and  Hypergeometric  Distribu¬ 
tion. 

Definition: 


Qn{x]  a,  f3,  N)  =  ^F2{—n,n  +  a  +  P+l,—x]a+l,—N]l),  n  =  0,l,...,N.  (A. 30) 
Orthogonality:  For  a  >  — 1  and  /?  >  — 1  or  for  a  <  —N  and  (3  <  —N, 

^ /a  +  a;\  //3  +  A  ^  (A.31) 

x=0  V  ^  /  \  ^  / 

where 

2  _  (~1)"(’2  +  Oi  +  (3  +  l)Ar+l(/3  +  l)n’2! 

"  (2n  +  a  +  /3  +  1)(q:  +  l)n( — A)„A! 

Recurrence  relation: 

—  xQn{x)  =  AnQn+lix)  —  (A„  +  C„)Q„(a;)  +  CnQn-lix),  (A. 32) 

where 


Qnix)  . —  Qn^Xj  O,  /?,  A) 


THE  WIENER-ASKEY  POLYNOMIAL  CHAOS 


25 


and 


_  {n+a+l3+l)(n+a+l){N—n) 

“  (2n+a+/3+l)(2n+a+/3+2) 


_  n(n+a+/3+Af+l)(n+/3) 
“  (2n+a+/3)(2n+a+/3+l)  ■ 


Rodriguez  formula: 


w(x;  a,  l3,  N)Qn{x]  a,  l3,  N) 


(-ir(/3+i). 

{-N)n 


V"[w(a;;  a  +  n,/3  +  n,  —  n)], 


(A.33) 


where 


w{x;  a,  /?,  N)  = 


a  +  a;^  f  P  +  N  —  x 

X 


N-x 

If  we  set  a  =  —a  —  1  and  /?  =  — /3  —  1,  we  obtain 

1  (JX/l.) 


= 


Apart  from  the  constant  factor  l/(  ^ 

distribution. 


^  ^  ,  this  is  the  definition  of  hyper  geometric 
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